# Paper: Is it still the economy? Economic voting in polarized politics
# County-level model
# Replication: Figure S16 (Supplementary Material)
remove(list = ls())

library(conflicted)
library(tidyverse)
conflict_prefer("filter","dplyr")
library(ggplot2)
library(dplyr)
library(here)
conflict_prefer("here", "here")
library(systemfit)
library(sandwich)
library(fixest)

load(here("Final_Paper", "R&R", "CountyLevel.RData"))



summary(ols_model <- lm(inc_share ~ unemp_rate*elitepolar_house + 
                          log(defl_pcincome) +
                          share_black + share_hisp + share_asian + 
                          share_young + share_elder + log(tot_pop) + 
                          log1p(share_college) + log1p(share_manuf) +
                          as.factor(rural_urban) + state_partiesdiff + 
                          lag_incshare + state + rep_inc + incumbent_cand, 
                        data = model_data))

coeffs_cl <- coeftest(ols_model, vcov = vcovCL, cluster = ~ county)
coeffs_cl

(est2 <- cbind(Estimate = coef(coeffs_cl), confint(coeffs_cl, level = 0.95)))


length(coeffs_cl[,1])
coeffs_cl[2,]
coeffs_cl[74,]

olsunem <- coeffs_cl[2,] #change in unemployment rate
olsint <- coeffs_cl[74,] #interaction

variables <- c("Unemp.", "Unemp.*Polarization")
ols_results <- rbind(olsunem, olsint)
glimpse(ols_results)

ols_data <- data.frame(variables, ols_results) %>% 
  mutate(model = "OLS") %>% 
  mutate(equation = "Vote Share")
glimpse(ols_data)

coef(summary(ols_model))[2,2]
coef(summary(ols_model))[74,2]    
se_unemp <- coef(summary(ols_model))[2,2]
se_int <- coef(summary(ols_model))[74,2]
se_ols <- rbind(se_unemp, se_int)
ols_data <- cbind(ols_data, se_ols)
glimpse(ols_data)

vcov(ols_model)[2, 2]
sqrt(vcov(ols_model)[2, 2])
sqrt(vcov(ols_model)[74, 74])
cov_unemp <- vcov(ols_model)[2, 74]


ols_data <- ols_data %>% 
  mutate(covb = ifelse(variables == "Unemp.*Polarization", cov_unemp, 0))

ols_data05 <- ols_data %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.35*Estimate, Estimate)) %>% 
  mutate(polarization = 0.35) 

ols_data06 <- ols_data %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.6*Estimate, Estimate)) %>% 
  mutate(polarization = 0.6) 
ols_data08 <- ols_data %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.85*Estimate, Estimate)) %>% 
  mutate(polarization = 0.85)

ols_data <- ols_data %>% 
  mutate(polarization = 1)

ols_data <- rbind(ols_data, ols_data05, ols_data06, ols_data08)
glimpse(ols_data)

ols_data <- ols_data %>% 
  mutate(varb = ifelse(variables == "Unemp.",
                       se_ols^2, polarization^2*se_ols^2)) 

ols_data2 <- ols_data %>% 
  mutate(variables = ifelse(variables == "Unemp.", "Unemp.*Polarization", variables)) %>% 
  group_by(equation, polarization, variables) %>% 
  summarise(Estimate = sum(Estimate),
            vars = sum(varb),
            covb = sum(covb)) %>% 
  ungroup()
glimpse(ols_data2)

# standard error of marginal effects
ols_data2 <- ols_data2 %>% 
  mutate(se = sqrt(vars + 2*polarization*covb))
ols_data2 <- ols_data2 %>% 
  select(equation:Estimate, se)

ols_data <- ols_data %>% 
  filter(polarization == 1) %>% 
  filter(variables == "Unemp.") %>% 
  mutate(polarization = 0)
ols_data <- ols_data %>% 
  select(equation, variables, polarization, Estimate, se = se_ols)

ols_data <- rbind(ols_data, ols_data2)
glimpse(ols_data)

ols_data <- ols_data %>% 
  mutate(icii = Estimate+(1.96*se)) %>% 
  mutate(ici = Estimate-(1.96*se))

sd(model_data$unemp_rate)
mean(model_data$unemp_rate)

sd(model_data$unemp_rate) #2.67
ols_data <- ols_data %>% 
  mutate(mean_unemp = mean(model_data$unemp_rate)) %>% 
  mutate(sd_uenmp = sd(model_data$unemp_rate)) %>% 
  mutate(estimate_mean_unemp = sd_uenmp*Estimate) %>% 
  mutate(icii_new = estimate_mean_unemp+(1.96*se)) %>% 
  mutate(ici_new = estimate_mean_unemp-(1.96*se))
glimpse(ols_data)


ols_data <- ols_data %>%
  mutate(polarization = as.factor(polarization)) %>%  
  mutate(polarization = factor(polarization, 
                               levels = c("0",
                                          "0.35",
                                          "0.6",
                                          "0.85",
                                          "1"))) %>% 
  mutate(equation = "Incumbent Party")

fig_se <- ols_data %>% 
  filter(polarization == "0.6" | polarization == "0.85") %>% 
  #filter(variables == "Unemp.Rate" | variables == "Unemp.Rate*Polarization") %>% 
  ggplot(aes(color = polarization, shape = polarization, fill = polarization)) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, size = .75) +
  geom_pointrange(aes(x = equation, y = estimate_mean_unemp, ymin = ici_new,
                      ymax = icii_new),
                  lwd = 0.8, position = position_dodge(width = .5)) + 
  theme_bw() + 
  ggtitle("(b) Clustered Standard Errors at the County Level") +
  ylab("Change in Vote Share") + xlab(NULL) + #ylim(-0.1, 0.05) +
  scale_color_manual(name = "Polarization:",
                     values = c("#009999", "#990000"), 
                     labels = c("0.60", "0.85")) +
  scale_shape_manual(name = "Polarization:",
                     values = c(24, 25), 
                     labels = c("0.60", "0.85")) +
  scale_fill_manual(name = "Polarization:",
                    values = c("#009999", "#990000"), 
                    labels = c("0.60", "0.85")) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 14))
fig_se


############### ############### ############### ############### 
############### ############### ############### ############### 
############### ############### ############### ############### 
############### County Fixed Effects
############### ############### ############### ############### 
############### ############### ############### ############### 
############### ############### ############### ############### 

summary(fe_model <- feols(inc_share ~ unemp_rate*elitepolar_house + 
                            log(defl_pcincome) +
                            share_black + share_hisp + share_asian + 
                            share_young + share_elder + log(tot_pop) + 
                            log1p(share_college) + log1p(share_manuf) +
                            as.factor(rural_urban) + state_partiesdiff + 
                            lag_incshare + rep_inc + incumbent_cand | county,
                          data = model_data))

(fe_inc <- cbind(Estimate = coef(fe_model), confint(fe_model, level = 0.95)))
fe_inc

fe_inc[1,]
fe_inc[25,]

summary(fe_model)
length(fe_model$coefficients)
fe_model$coefficients[1]
fe_model$coefficients[25]

olsunem <- fe_model$coefficients[1] #change in unemployment rate
olsint <- fe_model$coefficients[25] #interaction

olsunem <- fe_inc[1,] #change in unemployment rate
olsint <- fe_inc[25,] #interaction

variables <- c("Unemp.", "Unemp.*Polarization")
ols_results <- rbind(olsunem, olsint)
glimpse(ols_results)

fe_data <- data.frame(variables, ols_results) %>% 
  mutate(model = "FE") %>% 
  mutate(equation = "Vote Share")
glimpse(fe_data)

fe_model$se[1]
fe_model$se[25]  
se_unemp <- fe_model$se[1]
se_int <- fe_model$se[25]  
se_ols <- rbind(se_unemp, se_int)
se_ols
fe_data <- cbind(fe_data, se_ols)
glimpse(fe_data)

fe_data <- fe_data %>% 
  rename("se_ols" = "unemp_rate")

vcov(fe_model)
vcov(fe_model)[1, 1]
vcov(fe_model)[1, 25]
sqrt(vcov(fe_model)[1, 1])
sqrt(vcov(fe_model)[25, 25])
cov_unemp <- vcov(fe_model)[1, 25]


fe_data <- fe_data %>% 
  mutate(covb = ifelse(variables == "Unemp.*Polarization", cov_unemp, 0))

fe_data05 <- fe_data %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.35*Estimate, Estimate)) %>% 
  mutate(polarization = 0.35) 

fe_data06 <- fe_data %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.6*Estimate, Estimate)) %>% 
  mutate(polarization = 0.6) 
fe_data08 <- fe_data %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.85*Estimate, Estimate)) %>% 
  mutate(polarization = 0.85)

fe_data <- fe_data %>% 
  mutate(polarization = 1)

fe_data <- rbind(fe_data, fe_data05, fe_data06, fe_data08)
glimpse(fe_data)

fe_data <- fe_data %>% 
  mutate(varb = ifelse(variables == "Unemp.",
                       se_ols^2, polarization^2*se_ols^2)) 
glimpse(fe_data)

fe_data2 <- fe_data %>% 
  mutate(variables = ifelse(variables == "Unemp.", "Unemp.*Polarization", variables)) %>% 
  group_by(equation, polarization, variables) %>% 
  summarise(Estimate = sum(Estimate),
            vars = sum(varb),
            covb = sum(covb)) %>% 
  ungroup()
glimpse(fe_data2)

# standard error of marginal effects
fe_data2 <- fe_data2 %>% 
  mutate(se = sqrt(vars + 2*polarization*covb))
fe_data2 <- fe_data2 %>% 
  select(equation:Estimate, se)
glimpse(fe_data2)

glimpse(fe_data)
fe_data <- fe_data %>% 
  filter(polarization == 1) %>% 
  filter(variables == "Unemp.") %>% 
  mutate(polarization = 0)
fe_data <- fe_data %>% 
  select(equation, variables, polarization, Estimate, se = se_ols)

fe_data <- rbind(fe_data, fe_data2)
glimpse(fe_data)

fe_data <- fe_data %>% 
  mutate(se = ifelse(is.na(se), 0.001, se)) %>% 
  mutate(icii = Estimate+(1.96*se)) %>% 
  mutate(ici = Estimate-(1.96*se))

sd(model_data$unemp_rate)
mean(model_data$unemp_rate)

sd(model_data$unemp_rate) #2.67
fe_data <- fe_data %>% 
  mutate(mean_unemp = mean(model_data$unemp_rate)) %>% 
  mutate(sd_uenmp = sd(model_data$unemp_rate)) %>% 
  mutate(estimate_mean_unemp = sd_uenmp*Estimate) %>% 
  mutate(icii_new = estimate_mean_unemp+(1.96*se)) %>% 
  mutate(ici_new = estimate_mean_unemp-(1.96*se))
glimpse(fe_data)


fe_data <- fe_data %>%
  mutate(polarization = as.factor(polarization)) %>%  
  mutate(polarization = factor(polarization, 
                               levels = c("0",
                                          "0.35",
                                          "0.6",
                                          "0.85",
                                          "1"))) %>% 
  mutate(equation = "Incumbent Party")

fig_fe <- fe_data %>% 
  filter(polarization == "0.6" | polarization == "0.85") %>% 
  #filter(variables == "Unemp.Rate" | variables == "Unemp.Rate*Polarization") %>% 
  ggplot(aes(color = polarization, shape = polarization, fill = polarization)) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, size = .75) +
  geom_pointrange(aes(x = equation, y = estimate_mean_unemp, ymin = ici_new,
                      ymax = icii_new),
                  lwd = 0.8, position = position_dodge(width = .5)) + 
  theme_bw() + 
  ggtitle("(c) Model with County Fixed Effects") +
  ylab("Change in Vote Share") + xlab(NULL) + #ylim(-0.1, 0.05) +
  scale_color_manual(name = "Polarization:",
                     values = c("#009999", "#990000"), 
                     labels = c("0.60", "0.85")) +
  scale_shape_manual(name = "Polarization:",
                     values = c(24, 25), 
                     labels = c("0.60", "0.85")) +
  scale_fill_manual(name = "Polarization:",
                    values = c("#009999", "#990000"), 
                    labels = c("0.60", "0.85")) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 14))
fig_fe

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

############### ############### ############### ############### 
############### ############### ############### ############### 
############### ############### ############### ############### 


(est2 <- cbind(Estimate = coef(ols_model), confint(ols_model, level = 0.95)))


length(est2[,1])
est2[2,]
est2[74,]

olsunem <- est2[2,] #change in unemployment rate
olsint <- est2[74,] #interaction

variables <- c("Unemp.", "Unemp.*Polarization")
ols_results <- rbind(olsunem, olsint)
glimpse(ols_results)

ols_data2 <- data.frame(variables, ols_results) %>% 
  mutate(model = "OLS") %>% 
  mutate(equation = "Vote Share")
glimpse(ols_data2)

coef(summary(ols_model))[2,2]
coef(summary(ols_model))[74,2]    
se_unemp <- coef(summary(ols_model))[2,2]
se_int <- coef(summary(ols_model))[74,2]
se_ols <- rbind(se_unemp, se_int)
ols_data2 <- cbind(ols_data2, se_ols)
glimpse(ols_data2)

vcov(ols_model)[2, 2]
sqrt(vcov(ols_model)[2, 2])
sqrt(vcov(ols_model)[74, 74])
cov_unemp <- vcov(ols_model)[2, 74]


ols_data2 <- ols_data2 %>% 
  mutate(covb = ifelse(variables == "Unemp.*Polarization", cov_unemp, 0))

ols_data205 <- ols_data2 %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.35*Estimate, Estimate)) %>% 
  mutate(polarization = 0.35) 

ols_data206 <- ols_data2 %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.6*Estimate, Estimate)) %>% 
  mutate(polarization = 0.6) 
ols_data208 <- ols_data2 %>% 
  mutate(Estimate = ifelse(variables == "Unemp.*Polarization",
                           0.85*Estimate, Estimate)) %>% 
  mutate(polarization = 0.85)

ols_data2 <- ols_data2 %>% 
  mutate(polarization = 1)

ols_data2 <- rbind(ols_data2, ols_data205, ols_data206, ols_data208)
glimpse(ols_data2)

ols_data2 <- ols_data2 %>% 
  mutate(varb = ifelse(variables == "Unemp.",
                       se_ols^2, polarization^2*se_ols^2)) 

ols_data22 <- ols_data2 %>% 
  mutate(variables = ifelse(variables == "Unemp.", "Unemp.*Polarization", variables)) %>% 
  group_by(equation, polarization, variables) %>% 
  summarise(Estimate = sum(Estimate),
            vars = sum(varb),
            covb = sum(covb)) %>% 
  ungroup()
glimpse(ols_data22)

# standard error of marginal effects
ols_data22 <- ols_data22 %>% 
  mutate(se = sqrt(vars + 2*polarization*covb))
ols_data22 <- ols_data22 %>% 
  select(equation:Estimate, se)

ols_data2 <- ols_data2 %>% 
  filter(polarization == 1) %>% 
  filter(variables == "Unemp.") %>% 
  mutate(polarization = 0)
ols_data2 <- ols_data2 %>% 
  select(equation, variables, polarization, Estimate, se = se_ols)

ols_data2 <- rbind(ols_data2, ols_data22)
glimpse(ols_data2)

ols_data2 <- ols_data2 %>% 
  mutate(icii = Estimate+(1.96*se)) %>% 
  mutate(ici = Estimate-(1.96*se))

sd(model_data$unemp_rate)
mean(model_data$unemp_rate)

sd(model_data$unemp_rate) #2.67
ols_data2 <- ols_data2 %>% 
  mutate(mean_unemp = mean(model_data$unemp_rate)) %>% 
  mutate(sd_uenmp = sd(model_data$unemp_rate)) %>% 
  mutate(estimate_mean_unemp = sd_uenmp*Estimate) %>% 
  mutate(icii_new = estimate_mean_unemp+(1.96*se)) %>% 
  mutate(ici_new = estimate_mean_unemp-(1.96*se))
glimpse(ols_data2)


ols_data2 <- ols_data2 %>%
  mutate(polarization = as.factor(polarization)) %>%  
  mutate(polarization = factor(polarization, 
                               levels = c("0",
                                          "0.35",
                                          "0.6",
                                          "0.85",
                                          "1"))) %>% 
  mutate(equation = "Incumbent Party")

fig_olsa <- ols_data2 %>% 
  filter(polarization == "0.6" | polarization == "0.85") %>% 
  #filter(variables == "Unemp.Rate" | variables == "Unemp.Rate*Polarization") %>% 
  ggplot(aes(color = polarization, shape = polarization, fill = polarization)) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, size = .75) +
  geom_pointrange(aes(x = equation, y = estimate_mean_unemp, ymin = ici_new,
                      ymax = icii_new),
                  lwd = 0.8, position = position_dodge(width = .5)) + 
  theme_bw() + 
  ggtitle("(a) Unemployment Rate") +
  ylab("Change in Vote Share") + xlab(NULL) + #ylim(-0.1, 0.05) +
  scale_color_manual(name = "Polarization:",
                     values = c("#009999", "#990000"), 
                     labels = c("0.60", "0.85")) +
  scale_shape_manual(name = "Polarization:",
                     values = c(24, 25), 
                     labels = c("0.60", "0.85")) +
  scale_fill_manual(name = "Polarization:",
                    values = c("#009999", "#990000"), 
                    labels = c("0.60", "0.85")) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 14))
fig_olsa


multiplot(fig_olsa, fig_fe, fig_se, cols = 2)
